##################################################
## Author: Stephanie L. DeMora, Christian Lindke, Jenn Merolla, and Laura Stephenson
## Project: "Ready for a Woman President?" Paper Replication Script
## Date: September 25th 2020
##################################################

# Set up: Load in all packages. Please note that if "margins_summary" or "list.load" do not work, you may need to manually install "rlist" and "margins"
Packages <- c("ggplot2","hrbrthemes","tidyverse","haven","broom","Rmisc","dotwhisker","margins","jtools","mediation","here","rlist","dygraphs", "webshot", "htmlwidgets", "xts", "sensemakr")
lapply(Packages, library, character.only = TRUE)

# Plot Themes:
steph_theme <- list(
  hrbrthemes::theme_ipsum_rc(base_size = 12,grid= "",plot_title_size = 20,axis_text_size = 18),
  #ggsci::scale_color_material("grey"),
  #ggsci::scale_fill_material("grey"),
  theme(legend.position = "top", plot.caption = element_text(hjust = 0, face= "italic"),
        plot.title.position = "plot",
        plot.caption.position =  "plot"))

# Read in Data:
df <- read_dta("Data Files/CLEAN-CCES19_UCR.dta")

# Prep Data:
source("00. CCES_cleaning.R")


# ------------------------------------------------------------------------------------------------------------------------
# ----------------------------------------------  FIGURE 1  --------------------------------------------------------------
# ------------------------------------------------------------------------------------------------------------------------

# Harris:
var_a <- df1$UCR119b_1_1[df1$TREAT=="Control"]
var_b <- df1$UCR119b_1_1[df1$TREAT=="Neighbor"]
var_c <- df1$UCR119b_1_1[df1$TREAT=="Personal"]
var_d <- df1$UCR119b_1_1[df1$TREAT=="Harder Win"]
test1 <- t.test(var_b,var_a,conf.level = 0.90)
test2 <- t.test(var_c,var_a,conf.level = 0.90)
test3 <- t.test(var_d,var_a,conf.level = 0.90)

options(scipen = 999)
data <- tidy(lm(UCR119b_1_1 ~ TREAT, data=df1), conf.level = 0.90, conf.int = TRUE)
data$term <- c("(Intercept)","Neighbor","Personal","Harder Win")

tgc <- summarySE(df1, measurevar="UCR119b_1_1", groupvars=c("TREAT"), conf.interval = .90)
tgc$UCR119b_1_1 <- data$estimate
tgc$term <- tgc$TREAT
tgc$estimate <- tgc$UCR119b_1_1
tgc <- tgc[2:4,]

# two-tailed tests for all
tgc$conf.high <- c(test3[["conf.int"]][2], test1[["conf.int"]][2], test2[["conf.int"]][2])
tgc$conf.low <- c(test3[["conf.int"]][1], test1[["conf.int"]][1], test2[["conf.int"]][1])
tgc$conf.high_90 <- c(test3[["conf.int"]][2], test1[["conf.int"]][2], test2[["conf.int"]][2])
tgc$conf.low_90 <- c(test3[["conf.int"]][1], test1[["conf.int"]][1], test2[["conf.int"]][1])
test1 <- t.test(var_b,var_a,conf.level = 0.95)
test2 <- t.test(var_c,var_a,conf.level = 0.95)
test3 <- t.test(var_d,var_a,conf.level = 0.95)
tgc$conf.low_95 <- c(test3[["conf.int"]][1], test1[["conf.int"]][1], test2[["conf.int"]][1])
tgc$conf.high_95 <- c(test3[["conf.int"]][2], test1[["conf.int"]][2], test2[["conf.int"]][2])
tgc$model <- "Harris"
tgc1 <- tgc

# Warren:
var_a <- df1$UCR119d_1_1[df1$TREAT=="Control"]
var_b <- df1$UCR119d_1_1[df1$TREAT=="Neighbor"]
var_c <- df1$UCR119d_1_1[df1$TREAT=="Personal"]
var_d <- df1$UCR119d_1_1[df1$TREAT=="Harder Win"]
test1 <- t.test(var_b,var_a,conf.level = 0.90)
test2 <- t.test(var_c,var_a,conf.level = 0.90)
test3 <- t.test(var_d,var_a,conf.level = 0.90)
data <- tidy(lm(UCR119d_1_1 ~ TREAT, data=df1), conf.level = 0.90, conf.int = TRUE)
data$term <- c("(Intercept)","Neighbor","Personal","Harder Win")
tgc <- summarySE(df1, measurevar="UCR119d_1_1", groupvars=c("TREAT"), conf.interval = .90)
tgc$UCR119d_1_1 <- data$estimate
tgc$term <- tgc$TREAT
tgc$estimate <- tgc$UCR119d_1_1
tgc <- tgc[2:4,]
tgc$conf.low <- c(test3[["conf.int"]][1], test1[["conf.int"]][1], test2[["conf.int"]][1])
tgc$conf.high <- c(test3[["conf.int"]][2], test1[["conf.int"]][2], test2[["conf.int"]][2])
tgc$conf.low_90 <- c(test3[["conf.int"]][1], test1[["conf.int"]][1], test2[["conf.int"]][1])
tgc$conf.high_90 <- c(test3[["conf.int"]][2], test1[["conf.int"]][2], test2[["conf.int"]][2])
test1 <- t.test(var_b,var_a,conf.level = 0.95)
test2 <- t.test(var_c,var_a,conf.level = 0.95)
test3 <- t.test(var_d,var_a,conf.level = 0.95)
tgc$conf.low_95 <- c(test3[["conf.int"]][1], test1[["conf.int"]][1], test2[["conf.int"]][1])
tgc$conf.high_95 <- c(test3[["conf.int"]][2], test1[["conf.int"]][2], test2[["conf.int"]][2])
tgc$model <- "Warren"

tgc %>%
  dplyr::rename(UCR119b_1_1 = UCR119d_1_1) -> tgc
master <- rbind(tgc,tgc1)
master$TREAT <- car::recode(master$TREAT, "'Harder Win' = 'High Comfort 2'; 'Neighbor' = 'Low Comfort'; 'Personal' = 'High Comfort 1'")
master$term <- car::recode(master$term, "'Harder Win' = 'High Comfort 2'; 'Neighbor' = 'Low Comfort'; 'Personal' = 'High Comfort 1'")
master$TREAT <- fct_relevel(master$TREAT, "High Comfort 2","High Comfort 1","Low Comfort")

# Plot it:
ggplot(master, aes(x=TREAT, y=UCR119b_1_1, group=model, color=model)) +  
  geom_point(aes(x = TREAT, 
                 y = UCR119b_1_1), position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(x = TREAT, 
                     ymin = conf.low_90,
                     ymax = conf.high_90), position = position_dodge(width = 0.5),
                 lwd = 2) +
  geom_linerange(aes(x = TREAT, 
                     ymin = conf.low_95,
                     ymax = conf.high_95), position = position_dodge(width = 0.5),
                 lwd = 1/2) +
  scale_color_grey(start = 0.8, end = 0.2) + # Remove to bring color back
  labs(title = "Democratic Candidate Perceived Electability",
       subtitle = "Democratic Women Vs. Trump",
       y = "90% Confidence Intervals bolded, 95% Confidence Intervals unbolded",
       x = "Treatment",
       fill = "",
       caption = "Fall 2019 CCES Survey",
       color = "Candidate") +
  coord_flip() +
  steph_theme + 
  geom_hline(yintercept = 0, colour = "grey60", linetype = 2) -> Figure_1

Figure_1

# ------------------------------------------------------------------------------------------------------------------------
# ----------------------------------------------  FIGURE 2  --------------------------------------------------------------
# ------------------------------------------------------------------------------------------------------------------------

# Figure 2 Panel A ----------

# This is the code to produce our results from scratch. However, the results may differ slightly each time due to differing simulations.
# For this reason, we have included the .Rdata file to reproduce our identical results as well.

#set.seed(12082020)
# df2 <- df[!is.na(df$Gender_Vote),]
# model.0 <- lm(Gender_Vote ~ TREAT, df2)
# model.M <- lm(UCR119b_1_1 ~ TREAT, df2)
# model.Y <- lm(Gender_Vote ~ TREAT + UCR119b_1_1, df2)
# results4 <- mediate(model.M, model.Y, treat='TREAT', mediator='UCR119b_1_1',
#                     boot=TRUE, sims=1000, conf.level = .90, treat.value = "Neighbor", control.value = "Control")
#list.save(results4, file="Figure4_Results.RData")

results4 <- list.load("Data Files/Figure4_Results.RData")
plot(results4,sub = "90% Confidence Intervals", cex.sub = 1, cex.axis = 1.27, adj = .5, main = "Panel A: Harris")
summary(results4)

# Figure 2 Panel B ------------

# This is the code to produce our results from scratch. However, the results will differ each time due to differing simulations.
# For this reason, we have included the Rdata file to reproduce our identical results as well.

# set.seed(12082020)
# model.0 <- lm(Gender_Vote ~ TREAT, df2)
# model.M <- lm(UCR119d_1_1 ~ TREAT, df2)
# model.Y <- lm(Gender_Vote ~ TREAT + UCR119d_1_1, df2)
# results5 <- mediate(model.M, model.Y, treat='TREAT', mediator='UCR119d_1_1',
#                    boot=TRUE, sims=1000, conf.level = .90, treat.value = "Neighbor", control.value = "Control")
#list.save(results5, file="Figure5_Results.RData")

results5 <- list.load("Data Files/Figure5_Results.RData")
plot(results5,sub = "90% Confidence Intervals", cex.sub = 1, cex.axis = 1.27, adj = .5,main = "Panel B: Warren")
summary(results5)

# ------------------------------------------------------------------------------------------------------------------------
# ----------------------------------------------  FIGURE 3  --------------------------------------------------------------
# ------------------------------------------------------------------------------------------------------------------------

df1$party <- as.numeric(as.character(labelled::remove_var_label(df1$party)))

mod1 <- lm(UCR119b_1_1 ~ TREAT * GENDER * party, data = df1)
summary(mod1)
# calculate margins at 90% instead of 95%
margins_summary <- function(model, ..., level = 0.90, by_factor = TRUE) {
  mar <- margins(model, ...)
  summary(mar, level = 0.90, by_factor = by_factor)
}
tgc <- margins_summary(mod1, at = list(GENDER = range(df1$GENDER), party = range(df1$party)))
tgc$GENDER <- car::recode(tgc$GENDER, "0 = 'Men'; 1 = 'Women'")
tgc$party <- as.character(tgc$party)
tgc$model <- car::recode(tgc$party, "0 = 'Independent'; 1 = 'Democrat'")
tgc$term <- tgc$factor
tgc$conf.low = tgc$lower
tgc$conf.high = tgc$upper
tgc$conf.low_90 = tgc$lower
tgc$conf.high_90 = tgc$upper
tgc$estimate <- tgc$AME

# calculate margins at 95%
margins_summary2 <- function(model, ..., level = 0.95, by_factor = TRUE) {
  mar <- margins(model, ...)
  summary(mar, level = 0.95, by_factor = by_factor)
}

tgc95 <- margins_summary2(mod1, at = list(GENDER = range(df1$GENDER), party = range(df1$party)))

tgc$conf.high_95 <- tgc95$upper
tgc$conf.low_95 <- tgc95$lower
tgc1 <- subset(tgc,tgc$GENDER == "Women")
tgc1$factor <- paste(tgc1$factor, " ", sep = "")
tgc1$term <- paste(tgc1$term, " ", sep = "")
tgc <- subset(tgc,tgc$GENDER == "Men")

master <- rbind(tgc,tgc1)
master %>%
  filter(str_detect(term, "TREAT")) -> master

master$term <- substring(master$term, 6)
master %>%
  arrange(GENDER) -> master

master$term <- car::recode(master$term, "'Harder Win' = 'High Comfort 2'; 'Neighbor' = 'Low Comfort'; 'Personal' = 'High Comfort 1'")
master$term <- car::recode(master$term, "'Harder Win ' = 'High Comfort 2 '; 'Neighbor ' = 'Low Comfort '; 'Personal ' = 'High Comfort 1 '")

master$sort <- ""
master$sort[master$term == "Low Comfort"] <- 1
master$sort[master$term == "Low Comfort "] <- 4
master$sort[master$term == "High Comfort 1"] <- 2
master$sort[master$term == "High Comfort 1 "] <- 5
master$sort[master$term == "High Comfort 2"] <- 3
master$sort[master$term == "High Comfort 2 "] <- 6

master %>%
  arrange(sort) -> master

master$term <- factor(master$term,levels = c("High Comfort 2 ","High Comfort 1 ","Low Comfort ","High Comfort 2","High Comfort 1","Low Comfort"))

# Plot it:
ggplot(master, aes(x = term, y=AME, group=model, color=model)) +
  geom_point(aes(x = term, 
                 y = AME), position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(ymin = conf.low_90,
                     ymax = conf.high_90), position = position_dodge(width = 0.5),
                 lwd = 2) +
  geom_linerange(aes(ymin = conf.low_95,
                     ymax = conf.high_95), position = position_dodge(width = 0.5),
                 lwd = 1/2) +
  ggpubr::geom_bracket(data = master, xmin = 0.7, xmax = 3.5, y.position = 17,label = "Women",coord.flip = T,family = "Roboto") +
  ggpubr::geom_bracket(data = master, xmin = 3.7, xmax = 6.5, y.position = 17,label = "Men",coord.flip = T,family = "Roboto") +
  scale_color_grey(start = 0.8, end = 0.2) + # Remove to bring color back
  labs(title = "Democratic Candidate Perceived Electability",
       subtitle = "Harris Vs. Trump",
       caption = "Fall 2019 CCES Survey",
       x = "Treatment",
       y = "90% Confidence Intervals bolded, 95% Confidence Intervals unbolded",
       color = "Party",
       fill = "") +
  coord_flip() +
  steph_theme + geom_hline(yintercept = 0, colour = "grey60", linetype = 2) -> Figure_3

Figure_3

# ------------------------------------------------------------------------------------------------------------------------
# ----------------------------------------------  FIGURE 4  --------------------------------------------------------------
# ------------------------------------------------------------------------------------------------------------------------

mod1 <- lm(UCR119d_1_1 ~ TREAT * GENDER * party, data = df1)
summary(mod1)
huxtable::huxreg(mod1)

tgc <- margins_summary(mod1, at = list(GENDER = range(df1$GENDER), party = range(df1$party)))
tgc$GENDER <- car::recode(tgc$GENDER, "0 = 'Men'; 1 = 'Women'")
tgc$party <- as.character(tgc$party)
tgc$model <- car::recode(tgc$party, "0 = 'Independent'; 1 = 'Democrat'")
tgc$term <- tgc$factor
tgc$conf.low = tgc$lower
tgc$conf.high = tgc$upper
tgc$conf.low_90 = tgc$lower
tgc$conf.high_90 = tgc$upper
tgc$estimate <- tgc$AME

tgc95 <- margins_summary2(mod1, at = list(GENDER = range(df1$GENDER), party = range(df1$party)))
tgc$conf.high_95 <- tgc95$upper
tgc$conf.low_95 <- tgc95$lower
tgc1 <- subset(tgc,tgc$GENDER == "Women")
tgc1$factor <- paste(tgc1$factor, " ", sep = "")
tgc1$term <- paste(tgc1$term, " ", sep = "")
tgc <- subset(tgc,tgc$GENDER == "Men")
master <- rbind(tgc,tgc1)

master %>%
  filter(str_detect(term, "TREAT")) -> master

master$term <- substring(master$term, 6)

master %>%
  arrange(GENDER) -> master


master$term <- car::recode(master$term, "'Harder Win' = 'High Comfort 2'; 'Neighbor' = 'Low Comfort'; 'Personal' = 'High Comfort 1'")
master$term <- car::recode(master$term, "'Harder Win ' = 'High Comfort 2 '; 'Neighbor ' = 'Low Comfort '; 'Personal ' = 'High Comfort 1 '")


master$sort <- ""
master$sort[master$term == "Low Comfort"] <- 1
master$sort[master$term == "Low Comfort "] <- 4
master$sort[master$term == "High Comfort 1"] <- 2
master$sort[master$term == "High Comfort 1 "] <- 5
master$sort[master$term == "High Comfort 2"] <- 3
master$sort[master$term == "High Comfort 2 "] <- 6

master %>%
  arrange(sort) -> master
master$term <- factor(master$term,levels = c("High Comfort 2 ","High Comfort 1 ","Low Comfort ","High Comfort 2","High Comfort 1","Low Comfort"))

# Plot it:
ggplot(master, aes(x = term, y=AME, group=model, color=model)) +
  geom_point(aes(x = term, 
                 y = AME), position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(ymin = conf.low_90,
                     ymax = conf.high_90), position = position_dodge(width = 0.5),
                 lwd = 2) +
  geom_linerange(aes(ymin = conf.low_95,
                     ymax = conf.high_95), position = position_dodge(width = 0.5),
                 lwd = 1/2) +
  ggpubr::geom_bracket(data = master, xmin = 0.7, xmax = 3.5, y.position = 17,label = "Women",coord.flip = T,family = "Roboto") +
  ggpubr::geom_bracket(data = master, xmin = 3.7, xmax = 6.5, y.position = 17,label = "Men",coord.flip = T,family = "Roboto") +
  scale_color_grey(start = 0.8, end = 0.2) + # Remove to bring color back
  labs(title = "Democratic Candidate Perceived Electability",
       subtitle = "Warren Vs. Trump",
       caption = "Fall 2019 CCES Survey",
       x = "Treatment",
       y = "90% Confidence Intervals bolded, 95% Confidence Intervals unbolded",
       color = "Party",
       fill = "") +
  coord_flip() +
  steph_theme + geom_hline(yintercept = 0, colour = "grey60", linetype = 2) -> Figure_4

Figure_4

# ------------------------------------------------------------------------------------------------------------------------
# ----------------------------------------------  FIGURE 5  --------------------------------------------------------------
# ------------------------------------------------------------------------------------------------------------------------


linked <- lm(UCR119a_1_1 ~ TREAT, data=df1)
summary(linked)

data <- tidy(linked, conf.level = 0.90, conf.int = TRUE)
data$term <- c("(Intercept)","Neighbor","Personal","Harder Win")
data$model <- "Biden"
data2 <- tidy(linked, conf.level = 0.95, conf.int = TRUE)
data$conf.high_95 <- data2$conf.high
data$conf.low_95 <- data2$conf.low
data$conf.high_90 <- data$conf.high
data$conf.low_90 <- data$conf.low
tgc1 <- data
linked <- lm(UCR119c_1_1 ~ TREAT, data=df1)
summary(linked)
data <- tidy(linked, conf.level = 0.90, conf.int = TRUE)
data$term <- c("(Intercept)","Neighbor","Personal","Harder Win")
data$model <- "Sanders"
data2 <- tidy(linked, conf.level = 0.95, conf.int = TRUE)
data$conf.high_95 <- data2$conf.high
data$conf.low_95 <- data2$conf.low
data$conf.high_90 <- data$conf.high
data$conf.low_90 <- data$conf.low
tgc <- data
tgc <- tgc[2:4,]
tgc1 <- tgc1[2:4,]
master <- rbind(tgc,tgc1)
master$term <- car::recode(master$term, "'Harder Win' = 'High Comfort 2'; 'Neighbor' = 'Low Comfort'; 'Personal' = 'High Comfort 1'")
master$term <- factor(master$term,levels = c("High Comfort 2 ","High Comfort 1 ","Low Comfort ","High Comfort 2","High Comfort 1","Low Comfort"))

# Plot it:
ggplot(master, aes(x=term, y=estimate, group=model, color=model)) +
  geom_point(aes(x = term, 
                 y = estimate), position = position_dodge(width = 0.5)) + 
  geom_linerange(aes(x = term, 
                     ymin = conf.low_90,
                     ymax = conf.high_90), position = position_dodge(width = 0.5),
                 lwd = 2) +
  geom_linerange(aes(x = term, 
                     ymin = conf.low_95,
                     ymax = conf.high_95), position = position_dodge(width = 0.5),
                 lwd = 1/2) +
  scale_color_grey(start = 0.8, end = 0.2) + # Remove to bring color back
  labs(title = "Democratic Candidate Perceived Electability",
       subtitle = "Democratic Men Vs. Trump",
       y = "90% Confidence Intervals bolded, 95% Confidence Intervals unbolded",
       x = "Treatment",
       fill = "",
       caption = "Fall 2019 CCES Survey",
       color = "Candidate") +
  coord_flip() +
  steph_theme + geom_hline(yintercept = 0, colour = "grey60", linetype = 2) -> Figure_5

Figure_5
